Magnetism of Finite Graphene Samples: Mean-Field Theory compared with Exact 
Diagonalization and Quantum Monte Carlo Simulation 
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The magnetic properties of graphene on finite geometries are studied using a self-consistent mean- 
field theory of the Hubbard model. This approach is known to predict ferromagnetic edge states 
close to the zig-zag edges in single-layer graphene quantum dots and nanoribbons. In order to assess 
the accuracy of this method, we perform complementary exact diagonalization and quantum Monte 
Carlo simulations. We observe good quantitative agreement for all quantities investigated provided 
that the Coulomb interaction is not too strong. 

PACS numbers: 71.10.Fd; 73.22.Pr; 75.40.Mg 
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I. INTRODUCTION 

Graphene consists of a single layer of carbon atoms ar- 
ranged in a honeycomb crystal lattice^ and is a promising 
material with unique electronic properties. Among the 
most important characteristics, one should mention the 
presence of massless carriers, weak spin-orbit coupling, 
insensitivity to an external electrostatic potential (Klein 
paradox), fractional quantum Hall effect, etc. (for a re- 
view of the main features of graphene see Ref. [3)- The 
electronic properties of graphene nanostructures such as 
nanoribbons or quantum dots are expected to be very dif- 
ferent from bulk graphene. In fact, the Coulomb interac- 
tion is considerably enhanced in smaller geometries such 
as quantum dots, leading for example to unusual blockade 
effects.—"— On the other hand, the edge effect, which de- 
pends strongly on the geometry of the sample boundary, 
modifies the electronic structure of graphene /S^-S. In par- 
ticular, it has been predicted that finite graphene sam ples 
can exhibit magnetic edge states (see, e.g., Refs. I9l-[l9l) 
suggesting potential spintronics applications of graphene 
nanodevicesi^ 

It is common practice to use a mean-field theory 
(MFT) of the Hubbard model to investigate the mag- 
netic properties of graphene in finite geometries (see, 
e.g., Refs.iliSIiM^- Such a MFT is apphcable to any 
interaction and any geometry in a quite economic way: 
within the self-consistent mean-field (MF) approximation 
the main numerical effort is to solve the single- electron 
problem on a finite lattice. However, as far as we are 
aware, very little is known about the accuracy of this 
approximation. 

The main purpose of the present paper is to address 
this issue and check the accuracy of the MFT. We start 
by recalling a real-space formulation of the MFT in 
Sec. ini In Sec. [Hi] we briefly look at periodic bound- 
ary conditions'^ and show that we can reproduce edge- 



ferromagnetism for a dot with zig-zag edgesj -'^'^d^'^^ The 
accuracy of the MFT is carefully examined in Sec. IIVI 
where we present a comparison with exact diagonaliza- 
tion for a small "dot" and quantum Monte Carlo (QMC) 
simulations on a larger system with periodic boundary 
conditions. We conclude with a summary and perspec- 
tives in Sec. El 



II. MODEL AND COMPUTATION 

Since we are interested in the magnetic properties of 
graphene, interactions should be taken into account. To 
this end, we study the Hubbard model whose Hamilto- 
nian reads 



with rij 



are nearest neighbors on a 
honeycomb lattice. We denote the total number of sites 
by TV and the number of electrons with a spin projection 
(T=t,;by iV^. 

Due to the exponential growth of the Hilbert space 
dimension with N, a direct exact diagonalization of the 
Hubbard model (12. ip at half filling is only possible for 
system sizes until about 20 sites. In order to deal with 
larger system sizes we use a MF approximation: 
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MF 



-t E (2-2) 



It should be noted that the MF approximation breaks the 
SU(2)-symmetry of the original Hubbard model (|2.ip . 
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We compute the ground state 

\GS)= n < n 4/^io)' rf-,"=EQ-.»c-.i 

(2.3) 

and the one-electron spectrum eo-.a of the MF Hamiho- 
nian (1^^ using the LAPACK hbrary. 

This yields the ground-state energy, the local density 



^ ^ Qa,aiQiy,ia 



(2.4) 



the local magnetization (5,f) = ^ ('^i,! ~ "-i.J,): well as 



the spin correlation functions 
4 ^ 



{SIS^) = t( X! X! X! *5TjaQ^,»/5 



^Q=l p=l 

+Q^,iaQi,j0Ql^l3iQ\^ao I (2-5) 



for i ^ j and 
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(n,t) + (n«i)-2(n,t)(n4) (2.6) 



for i — j. 

Self-consistency requires that the expectation values 
(ricr,i) entering (|2.2p are equal to the expression (|2.4I) de- 
rived from this Hamiltonian. We solve this condition it- 
eratively using suitable initial conditions with given num- 
bers of electrons No-. To overcome convergence problems, 
we use a thermal state compatible with the Fermi-Dirac 
distribution at a given temperature instead of the ground 
state for the first iterations32i In this case the average 
density is computed as 



where fio- is a set of single-particle states chosen ran- 
domly with probability n{ea.a) = 1/(1 + el'^''^"-""'^)!). 



III. RESULTS OF THE MEAN-FIELD 
APPROXIMATION 

A. System with periodic boundary conditions 

First we briefly discuss the MFT for the infinite sys- 
tem with periodic boundary conditions. If we assume 
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FIG. 1: (Color online) Mean-field result for the edge magne- 
tization of a hexagonal graphene quantum dot with A'' = 96 
sites and zig-zag edges (at half filling and with U — 2t). 



a Neel-ordered configuration, we find a Mott-Hubbard 
phase transition at the literature value Uc ~ 2.23 1,^^ 
where the system goes from a paramagnetic semi-metal 
to an insulator with antiferromagnetic order. The asymp- 
totic behavior of the Neel order parameter and the single- 
particle gap is numerically consistent with linear behavior 
m U — Uc for U > Uc, i.e., associated critical exponents 
equal to one. These unusual mean-field exponents refiect 
the unusual density of states^ which is linear close to 
the Fermi energy (compare also Ref. |23[ ) . 

We use the critical value Uc mainly to choose an ap- 
proximate value of U to describe graphene. In fact, the 
correct value of the Coulomb interaction in graphene is 
not yet known. Taking the value of U in polyacetylene, 
where U = lOeV and t = 2.5eV, suggests [/ « 4i for 
graphenei^ Since, at the MF level, this value locates the 
system well inside the antiferromagnetic phase and it is 
observed that large graphene sheets do not show mag- 
netic order, we have decided to use a value of U smaller 
than Uc, U = 2t for the following computation. 



B. Edge magnetism on zig-zag edge 

It is well known that even for values of U smaller than 
the critical value Uc, one observes a form of ferromag- 
netism on the zig-zag edge of a graphene ribbon^ii^ii^iii 
or a quantum dotJ^ii^iiS As an example, Fig. [T] shows 
our results for the local magnetization of a hexagonal dot 
with 96 sites. One observes local ferromagnetic behavior 
at each zig-zag edge. By contrast, systems with armchair 
edges do not show specific magnetic properties and follow 
an evolution closer to the one of a system with periodic 
boundary conditions. The difference between the two 
edges appears to be a consequence of the fact that in the 
zig-zag case only one sublattice is represented on the edge 
while in the armchair case both sublattices are present. 
Detailed explorationsi^iii demonstrated that the ferro- 
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FIG. 2: (Color online) Comparison MFT-ED for the finite- 
size system of 16 sites sketched in the inset of the top panel 
at half filling. The bottom panel shows the total staggered 
magnetization Ms in the main panel and the z-component Al^ 
in the inset. 



magnetism of zig-zag edges resists to armchair defects 
and appears already for short edges. 



IV. ACCURACY OF THE APPROXIMATION 

A. Comparison with exact diagonalization for open 
boundary conditions 

To verify the accuracy of the MFT wc first compare the 
results with those obtained with exact diagonalization 
(ED) of the Hubbard model which was performed with 
Spinpackj^ Due to the exponential growth of the Hilbert 
space, ED is limited to very small systems. Here we have 
studied the dot-like cluster of 16 sites shown in the inset 
of the top panel of Fig. [31 The following quantities were 
computed: (i) The ground-state energy, (ii) The charge 
gap defined as AE = E^-i — 2En + -Eat+i, where En 
is the ground-state energy in the sector with n electrons. 
Within MFT the charge gap can be identified with the 



single-particle gap Agp. (iii) The z- and total staggered 
magnetization as defined in terms of the longitudinal and 
total spin structure factor 



= (4.1) 
= ^ ^J^iy70~^) . (4.2) 



Here, (—1)*^-' is a short-hand notation for +1 (—1) if i 
and j belong to the same (different) sublattice(s). Within 
MFT, the correlation functions appearing in (I4.ip and 
()4.2p are computed from p.Sp and p.6|) . In a numerical 
solution of the Hubbard model respecting SU(2) symme- 
try one finds = Ms/V^- 

Figure [5] shows a comparison of ground-state energy, 
charge gap, and the two staggered magnetizations com- 
puted both with MFT and ED at half filling, i.e., a total 
of 16 electrons. As expected, the two methods yield iden- 
tical results ioi U — and the MF ground-state energy is 
always above the exact answer. The results for all three 
quantities stay close for U <2t. This supports the appli- 
cability of MFT at least as a semi-quantitative method in 
particular for the parameters of the dot shown in Fig. [T] 



B. Comparison to Quantum Monte Carlo for 
periodic boundary conditions 

In order to assess the quality of the MFT for larger 
but still finite systems, we employ QMC simulations. We 
use a projective determinantal QMC approach^ to ob- 
tain ground-state properties at half filling. Within this 
scheme, expectation values of a physical observable A are 
obtained from 



(A) = lim 

0— >CXD 



(^'T|e-e^|*T> 



(4.3) 



where the trial wavefunction |\1/t) must be non- 
orthogonal to the ground state and Q corresponds to a 
projection parameter. We found = 40/t to be suffi- 
cient to obtain converged ground-state quantities within 
the statistical uncertainty. For the presented simula- 
tions, 8 was split into discrete step At in the Trotter 
decomposition. We verified by extrapolating At — )■ 
that taking At — 0.05/t produced no discretization arti- 
facts. The simulations were performed on systems with 
periodic boundary conditions. 

Figure [3] shows our QMC and MFT data for a finite 
system with N — 162 lattice sites. The top panel of Fig. [3] 
shows the energy per site. In the middle panel of Fig. [3] 
we compare the single-particle gap Agp between MFT 
and QMC. In the QMC simulations Asp(fc) was obtained 
by fitting the exponential tail of the imaginary-time dis- 
placed Green's function G{k,T) oc exp(— TAsp(fc)) at 
large imaginary time t. The single-particle gap Agp 
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FIG. 3: (Color online) Comparison MFT-QMC for a system 
with periodic boundary conditions and A*' = 162 sites at half 
filling. QMC error bars are smaller than the size of the sym- 
bols. 



shown in Fig. |3] equals Asp(-ft'), i.e., the smallest excita- 
tion gap at the Dirac points. Finally, the bottom panel 
of Fig. [3] compares the QMC results for the the total 
staggered magnetization Eq. (|4.2p (main panel) and the 
z-component Eq. (|4.ip (inset) with the MFT results. 

Again, the MFT follows the QMC results closely for 
U <2t. In the present case, the MF curves for the en- 
ergy, single-particle gap Agp, and are always above 
the QMC curves. In fact, one observes that in the regime 
[/ < 2 1 the agreement between MFT and QMC is a bit 
better for than for M^- This can be attributed to the 
MF approximation explicitly breaking the SU(2) symme- 
try. Indeed, in this case one finds that the MF contribu- 
tion of the X- and y-components to Ms are independent 
of U for U > whereas the z-component increases with 
increasing U . 

Appreciable quantitative differences can be observed 
in Fig.[3]at large U in particular in Agp and Ms- Indeed, 
MFT is known to underestimate the stability-range of the 
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I.e., 



paramagnetic semi-metal by about a factor 2 
quantitative differences are expected for U larger than 



the mean-field critical value Uc « 2.23 i. Furthermore, 
in the limit [/ — >■ oo we expect to recover the S* = 1/2 
Heisenberg model where it is known (compare, e.g., Refs. 
'27','28) that a full quantum mechanical treatment of the 
quantity defined in Eq. (j4.2p yields a value which for 
large N is only about 55% of the classical (i.e., MF) 
value Ms = 1/2. The difference at the right boundary 
of the bottom panel of Fig. [3] is indeed of this order. 
While the MFT and QMC deviate in the precise posi- 
tion of the quantum critical point, they agree on locating 
the system in the paramagnetic semi-metallic phase for 
U < Uc ^ 2.23 1 and the correspondence is at least semi- 
quantitative for a finite-size system and U ^2t. 



V. CONCLUSION AND PERSPECTIVES 

We investigated a self-consistent mean-field approxi- 
mation to the Hubbard model on the honeycomb lattice, 
concentrating on half filling. The infinite system exhibits 
a Mott-Hubbard transition from a paramagnetic semi- 
metal for a Coulomb repulsion U < Uc to an antiferro- 
magnetic insulator for U > Uc with a MF critical value 
Uc ~ 2.23 The mean-field critical exponents associ- 
ated to the gap and Neel order parameter are numerically 
consistent with the value onei^i 

We studied the accuracy of the MFT for finite-size 
systems with complementary exact diagonalization and 
quantum Monte Carlo simulations of the Hubbard model. 
We computed the ground-state energy, the single-particle 
gap, and the staggered magnetizations obtained from the 
spin correlation functions. The MFT reproduces the 
qualitative behavior found by the other two methods. 
Furthermore, the quantitative agreement is reasonable 
for t/ < 2t, i.e., the region which is identified as a para- 
magnetic semi-metal both by MFT and QMC. For large 
values of U, quantitative differences become appreciable. 
However, the latter regime corresponds to an antiferro- 
magnetic insulator, not relevant to graphene. 

A weak-coupling instability to a canted antiferromag- 
net emerges in graphene when an in-plane magnetic field 
is turned oni^ This weak-coupling instability can be 
captured at the mean-field level and is confirmed by 
QMC simulations.^^ Being equally a weak-coupling phe- 
nomenon, we believe that it will be possible to observe 
edge ferromagnetism in future QMC simulations with zig- 
zag boundary conditions. Nevertheless, the MFT has ac- 
cess to bigger systems than QMC since it is numerically 
less demanding. Here we have presented the first piece 
of evidence that the MFT may be expected to be quan- 
titatively reliable for [/ < 2 1. 
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